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We construct quasiequilibrium sequences of black hole-neutron star binaries in general relativity. 
We solve Einstein's constraint equations in the conformal thin-sandwich formalism, subject to black 
hole boundary conditions imposed on the surface of an excised sphere, together with the relativistic 
equations of hydrostatic equilibrium. In contrast to our previous calculations we adopt a flat spatial 
background geometry and do not assume extreme mass ratios. We adopt a F = 2 polytropic 
equation of state and focus on irrotational neutron star configurations as well as approximately 
nonspinning black holes. We present numerical results for ratios of the black hole's irreducible mass 
to the neutron star's ADM mass in isolation of M^?j^/Madm,o = li 2, 3, 5, and 10. We consider 
neutron stars of baryon rest mass M^^/M§'^'^ = 83% and 56%, where is the maximum 

allowed rest mass of a spherical star in isolation for our equation of state. For these sequences, 
we locate the onset of tidal disruption and, in cases with sufficiently large mass ratios and neutron 
star compactions, the innermost stable circular orbit. We compare with previous results for black 
hole-neutron star binaries and find excellent agreement with third-order post-Newtonian results, 
especially for large binary separations. We also use our results to estimate the energy spectrum of 
the outgoing gravitational radiation emitted during the inspiral phase for these binaries. 
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I. INTRODUCTION 

Coalescing compact binaries composed of neutron stars 
and/or black holes are among the most promising sources 
of gravitational waves for ground-based laser interferom- 
eters [H [3, H, 13] ■ Observations of short gamma-ray bursts 
(SGRBs) by the Swift and HETE-2 satellites (see, e.g. Q 
and references therein) suggest that their central engine 
may well be the merger remnant of compact binaries that 
contain a neutron star, namely the remnants of black 
hole-neutron star (BHNS) binaries or binary neutron star 
(BNS) systems d 0, H ll]. Motivated by these consid- 
erations, significant theoretical effort has gone into the 
modeling of these binaries. 

So far, most studies of BHNS binaries have been per- 
formed within the framework of Newtonian gravity in ei- 
ther some or all aspects of the calculation (see, e.g. 
HH, [H, El, [3JllMJl3,J3 for quasiequilibrium cal- 
culations and |l9l. l2Cll. fell. l22l [23l. [24 . Hsj for dynamical 
simulations). More recently, several groups have initiated 
studies of BHNS binaries in a fully relativistic framework, 
both for quasiequilibrium models ^26", '2?!, [H, [2^, Is^] and 
dynamical simulations 7, 31, 32, 33, 34]. 

In our own effort, we initially studied BHNS quasiequi- 
librium models under the assumption of extreme mass 
ratios, where the mass of the black hole is much greater 
than that of the neutron star 0, [13, [H, [3l| . This assump- 
tion is very appealing computationally and therefore a 
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natural first step in a systematic study of these binaries, 
but breaks down for the astrophysically more interesting 
comparable-mass binaries. For the latter, the neutron 
star is tidally disrupted outside of the black hole's inner- 
most stable circular orbit (ISCO) (see, e.g. [Tsl]) and 
may form an accretion disk, which is considered crucial 
for launching of a gamma-ray burst. The observation of 
such a process would provide a wealth of astrophysical in- 
formation (35! . In a subsequent study of quasiequilibrium 
BHNS models ^29^], we therefore relaxed the assumption 
of extreme mass ratios and studied comparable-mass bi- 
naries (see also [sol. Is^]). 

We construct quasiequilibrium models by solving 
Einstein's constraint equations in the conformal thin- 
sandwich formalism (see, e.g. [H, [13, [HI ) . These equa- 
tions constrain only some of the gravitational fields; oth- 
ers, namely the conformal background solution, are freely 
specifiable and have to be chosen before the equations 
can be solved (see Section III Al below) . In our previous 
studies [23, [H, [2^ we chose a background solution de- 
scribing a single black hole in Kerr-Schild coordinates (in 
[28t we also include results for a spatially flat background 
for extreme mass ratio binaries). This seemed appealing 
in a number of ways; for example, it allows for black hole 
spin and the coordinate system extends smoothly into 
the black hole interior. However, in [2^ we found siz- 
able deviations from post-Newtonian (PN) results which 
we could attribute, at least in part, to the choice of the 
background solution [H, [s^ . Moreover, the recent black 
hole boundary conditions of Cook and Pfeiffer [s^ (see 
also [4(i]) allow for black hole rotation and penetrating 
coordinates even for a spatially flat background solution, 
which simplifies the equations significantly. 

In this paper we therefore revisit quasiequilibrium 
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models of BHNS binaries. We now use a spatially 
flat background solution instead of a Kerr-Schild metric 
(compare [sO]), and we also adopt a new decomposition 
of the field variables which significantly reduces the nu- 
merical error that was still present in [29| (see Appendix 
Ej) . We construct sequences of BHNS binaries in quasi- 
circular orbits for fixed black hole irreducible masses and 
neutron star baryon rest masses, and focus on irrota- 
tional neutron stars orbiting approximately nonspinning 
black holes (we will implement a more accurate condi- 
tion for nonspinning black holes in the future, following 
[40|). We find excellent agreement with post-Newtonian 
results for large binary separations. We track cusp forma- 
tion to locate the onset of tidal disruption, and, for large 
mass ratios and neutron star compactions, find turning- 
points on a binding energy curve to locate the ISCO. For 
those sequences that encounter an ISCO we find reason- 
able agreement with the ISCO of test particles around 
Schwarzschild black holes. 

The paper is organized as follows. We briefly review 
the basic equations in Section II. We present numerical 
results in Section III, and provide a comparison with the 
results of [H, (291, |30[ in Section IV. In Section V we 
briefly summarize our findings. 

Throughout this paper we adopt geometric units with 
G = c = 1, where G denotes the gravitational constant 
and c the speed of light. Latin and Greek indices denote 
purely spatial and spacetime components, respectively. 

II. FORMULATION 

In this Section we briefly review the equations we solve 
to obtain a BHNS binary in quasiequilibrium. For a more 
detailed discussion we refer to the review articles [13, [sll 
as well as Sec. II of [4l| for the hydrostatics. 

A. Gravitational field equations 

The line element in 3 + 1 form is written as 

= -a^dt^ + -fijidx' + (3'dt){dx^ + (3^dt), (1) 

where g^j/ is the spacetime metric, a the lapse function, 
/3' the shift vector, and 7^- the spatial metric induced 
on a spatial slice E. The spatial metric ^ij is further 
decomposed according to 7^ = ip^jij, where ip denotes 
the conformal factor and jij the background spatial met- 
ric. We also decompose the extrinsic curvature iC'^ into 
a trace K and a traceless part A^^ according to 

K'^ = ^/j-^°A'' + --f'^K. (2) 
3 

The Hamiltonian constraint then becomes 

= -27r^V -I- + ^tp^K^ - U~'A,jA'\ (3) 



Here = j^^^iVj is the covariant Laplace operator, 
Vi the covariant derivative, Rij the Ricci tensor, and 
R = 7*-'i?ij the scalar curvature, all associated with the 
conformal background metric "fij . 

We employ the conformal thin-sandwich decomposi- 
tion of the Einstein equations [IBl- In this decomposi- 
tion, we use the evolution equation for the spatial metric 
to express the traceless part of the extrinsic curvature in 
terms of the time derivative of the background metric, 
Uij = dtjij, and the gradients of the shift vector. Under 
the assumption of equilibrium, i.e., Uij = in a corotat- 
ing coordinate system, the traceless part of the extrinsic 
curvature reduces to 

^■U ^ f_ fY^pl + _ lf^y,(]A . (4) 

2a V 3 / 

Inserting Eq. ^ into the momentum constraint we ob- 
tain 

V^^'-f ^V*(V,/3^')+K/3^' 
3 

= leiraiP^f + 2A'^\/j{a'4j-^) + -af^\7jK. (5) 

3 

For the construction of quasiequilibrium data it is also 
reasonable to assume dtK = Q in a, corotating coordi- 
nate system. The trace of the evolution equation for the 
extrinsic curvature then yields 

V^a = A'Kai!'^{p + S) + \a^^K^ + %l^(i'\7iK 

+a^l)-^Ai^jA'^ - 2f^V,aVj In ^. (6) 

Equations ((3]), (O and ([6|) provide equations for the 
lapse function a, the shift vector and the conformal 
factor V', while A^^ can be found from Eq. (jlj. The con- 
formally related spatial metric 7^ and the trace of the 
extrinsic curvature K remain freely specifiable, and have 
to be chosen before we can solve the above equations 
(note that we have already set to zero the time deriva- 
tives of these quantities, which are also freely specifi- 
able). In [HI we identified these quantities with the cor- 
responding quantities for a Schwarzschild black hole in 
Kerr-Schild coordinates. Instead, we now assume a flat 
background 7^ = 77^, where rjij denotes a flat spatial 
metric, and maximal slicing K — 0. In Cartesian coordi- 
nates, Eqs. jS]) and ^ then reduce to 

Aij = -2n^''p - l^/j-'A,A'\ (7) 
o 

A/3' + \d\djl3^) ^ lenatP^f + 2A'^ djiai;'^), (8) 

Aa = Aira^j'^ip + S) + a^-^A,jA'^ 

--2r]'^diadj\n^, (9) 

where A and di denote the flat Laplace operator and the 
flat partial derivative, while Eq. ^ becomes 
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For numerical purposes we further decompose the vari- 
ables and their equations into parts associated with the 
black hole and the neutron star. For details of this de- 
composition we refer to Appendix |^ but we note that 
we have found this decomposition crucial for improving 
the accuracy of the calculations. 

The matter terms on the right-hand side of Eqs. (O, 
([5]), and (O are derived from the projections of the stress- 
energy tensor T^^ into the spatial slice E. Assuming an 
ideal fluid, we have 



T. 



(11) 



where is the fluid 4-velocity, po the baryon rest-mass 
density, pi the internal energy density, and P the pres- 
sure. Denoting the future-oriented unit normal to S as 



the relevant projections of T^i^ are 



(12) 
(13) 
(14) 
(15) 



B. Hydrostatic equations 



The matter in the neutron star interior has to satisfy 
the relativistic equations of hydrodynamics. For station- 
ary configurations, the relativistic Eulcr equation can be 
integrated once to yield 



7 

ha — = constant, 

70 



(16) 



where h = {po + Pi + P) / Po is the fluid speciflc enthalpy, 
and 7 and 70 are Lorenz factors between the fluid and the 
rotating frame, and the rotating frame and the inertial 
frame (see Sec. II. C. of [2^ for the definitions). For 
irrotational fluids, the fluid 3- velocity with respect to the 
inertial observer, can be expressed in terms of the 
gradient of a velocity potential ^ as 



(17) 



where u* is the time component of the fluid 4-velocity 
M^. Having taken into account the expression of the 3- 
velocity, the equation of continuity becomes 



(V^*)VJ^ 



Po 



0, 



(18) 



where denotes the covariant derivative associated 
with g^j,. We refer to ^1 for a more detailed deriva- 
tion of the hydrostatic equations, and to [H, IH, H^l for 
BNS apphcations. 



C. Equation of state 

We adopt a polytropic equation of state in the form 
P = npl, (19) 



where F denotes the adiabatic index and k a constant. 
In this paper, we focus on the case F = 2. Using this 
equation of state, the definition of the speciflc enthalpy 
h — [pq + pi + P)/p(), and a thermodynamic relation 
(Gibbs-Duhem relation) 



dh 
~h 



dP 



PQ+ Pi+ P 

we obtain the internal energy density 



F- 1 



(20) 



(21) 



Since dimensions enter the problem only through the 
constant k, it is convenient to rescale all dimensional 
quantities with respect to the polytropic length scale 



R 



poly 



a/(2r-2) 



(22) 



D. Boundary conditions 



In order to solve the gravitational fleld equations ([7]), 
(HI), and (HI), we have to set appropriate boundary con- 
ditions on two different boundaries: outer boundaries at 
spatial infinity and inner boundaries on the black hole 
horizons. 

The boundary conditions at spatial infinity follow from 
the assumption of asymptotic flatness. With the help of a 
radial coordinate transformation u = 1/r in the external 
computational domain our computational grid extends 
to spatial infinity 'Jl], |4^ , and we can impose the exact 
boundary conditions 



V'lr^oo = 1, 

/3%_oo ^{nx R)\ 



(23) 
(24) 
(25) 



where VL is the orbital angular velocity of the binary 
system measured at infinity, and R = [X, Y, Z) is a 
Cartesian coordinate which origin is located at the cen- 
ter of mass of the binary system. Here, we express the 
shift vector in a corotating coordinate system that we 
adopt throughout our calculation. In an inertial coordi- 
nate system, the shift vector would tend to zero at spatial 
infinity, while in the corotating coordinate system of the 
numerical code the shift vector diverges at spatial infin- 
ity. For computational purposes it is therefore convenient 
to write the shift vector as a sum of the rotational shift 
term /3*q^ = (17 x Ry and a residual part (which tends 
to zero at spatial inflnity), and solve the equations only 
for the latter (see Appendix 1X1 for a detailed explanation 
of the decomposition of the metric quantities and related 
equations). 

The inner boundary conditions arise from the excision 
of the black hole interior. The assumption that the black 
hole is in equilibrium leads to a set of boundary condi- 
tions for the conformal factor and shift vector 13911 (see 
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also [40, |46[ as well as the related isolated horizon for- 
malism, e.g. [U El]). The boundary condition for the 
conformal factor is 



(26) 



where s* = -0 P is the outward pointing unit vector 



normal to the excision surface and h 
metric on the excision surface, hij 



is the induced 



The quantity J is computed from the projection of the 
extrinsic curvature Kij as J = k'^Kij. The boundary 
condition on the normal component of the shift vector is 

P±\s^a\s, (27) 

while the tangential components have to satisfy 

Pi^ls ^ e],nix\ (28) 

Here fi;^ is the black hole spin angular velocity vector, 
which we take to be aligned with the Z-axis, and x'^ is a 
Cartesian coordinate centered on the black hole. To con- 
struct approximately nonspinning black holes, we set the 
shift vector according to Eqs. ([27| and (|28p with il.r — fl. 
This assignment corresponds to a "leading-order approx- 
imation" in the language of [40| . and we plan to imple- 
ment their more accurate condition for nonspinning black 
holes in the future. According to [1^ , the boundary con- 
dition on the lapse function can be chosen freely. In this 
paper, we choose a Neumann boundary condition 



da 
dr 







(29) 



on the excised surface, where r is a radial isotropic coor- 
dinate. 



E. Orbital angular velocity 

In this paper we take the rotational axis of the binary 
system to be the Z-axis, and the line connecting the black 
hole and neutron star centers to be the X-axis. To impose 
a quasicircular orbit we require a force balance along the 
X-axis at the center of the neutron star, which results in 
the condition 1411 



dlnh 



dX 



(XNS,yNS:0) 



= 0. 



(30) 



Here Xns and Yns denote the X and y-coordinates of 
the center of the neutron star relative to the rotational 
axis of the binary system. Imposing this condition de- 
termines the orbital angular velocity fl. We confirm 
that the orbital angular velocity obtained by this method 
agrees with that obtained by requiring the enthalpy at 
two points on the neutron star's surface to be equal [l^l 
within one part in 10^^. 



F. Global quantities 

It is reasonable to require that the irreducible mass of 
the black hole and the baryon rest mass of the neutron 
star are conserved during the inspiral of the BHNS bina- 
ries. For such a constant-mass sequence we then moni- 
tor the Arnowitt-Deser-Misner (ADM) mass, the Komar 
mass, and the total angular momentum, all of which are 
defined globally as follows. 

The irreducible mass of the black hole is defined as 



-"^^irr — 



A 



EH 



167r ' 



(31) 



where Aejj is the proper area of the event horizon. In 
practice we approximate this area with that of the appar- 
ent horizon, Aah , which is computed from an integral on 
the excision surface S, 



The baryon rest mass of the neutron star is 



NS 



pou'V^d^x, 



(32) 



(33) 



where g is the determinant of g^jy. In our case we have 
\/—g = a4>^, so that 



In the polytropic units of Eq. 
baryon rest mass according to 



NS _ ^'-'b 



poly 



(34) 

we normalize the 
(35) 



whereby Mg ^ is the rest mass of the polytrope for poly- 
tropic constant n = 1. 

We express the ADM mass in isotropic Cartesian co- 
ordinates as 



The Komar mass can be written as 

MKomar = ^ <^ 9'adS,, 



(36) 



(37) 



where we use the fact that the shift vector falls off suffi- 
ciently rapidly. The total angular momentum is 



1 

16^^ 



k i {X'K''^ ~ X^K^^)dSi, 



(38) 



where X* is a spatial Cartesian coordinate relative to the 
center of mass of the binary system. Finally, the linear 
momentum is 



(39) 
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where we have assumed maximal slicing K = 0. 

During the iteration, we require that the linear momen- 
tum vanishes. We enforce this by changing the position 
of the black hole and neutron star relative to the location 
of the axis of rotation. Stated differently, this condition 
fixes the location of the axis of rotation. 

We define the binding energy of the binary system as 



Ey, = Madm - Mo, 



(40) 



where Mq is the ADM mass of the binary system at in- 
finite orbital separation, as defined by the sum of the 
irreducible mass of the isolated black hole and the ADM 
mass of an isolated neutron star with the same baryon 
rest mass, 



Mo = M, 



BH 



■ -'"-'ADM.O- 



(41) 



In order to measure a global error in the numerical re- 
sults, we define the virial error as the fractional difference 
between the ADM mass and Komar mass, 

AfADM — -^^Komar 



5M 



Ma 



DM 



(42) 



Mg^ = 0.15 and 0.10, with corresponding compactions 
of Madm,o/^o = 0.145 and 0.0879. Here i?o is the areal 
radius of the neutron star in isolation. The maximum 
baryon rest mass for spherical F = 2 polytropes in isola- 
tion is M^^^ — 0.180. Our more compact polytrope with 
Mg^ = 0.15 has a compaction that is reasonably realis- 
tic (e.g. it would apply to binaries containing neutron 
stars with baryon rest mass 1.5A^0 and isolated spheri- 
cal neutron stars with a maximum baryon rest mass of 
I.8M0); we consider the less compact model for purposes 
of comparison. 

We also consider mass ratios M^^/M^y)M — 1' 2, 3, 
5, and 10. Note again that we fix the irreducible mass of 
the black hole and the baryon rest mass of the neutron 
star for the construction of constant-mass sequences. For 
the definition of the mass ratio, however, we use the ADM 
mass of a spherical isolated neutron star M^^j^ g, since 
this turns out to be more convenient for comparisons with 
third-order post- Newtonian results. 

We summarize our results in Tables IIVI and |V] in Ap- 
pendix [B] 



III. NUMERICAL RESULTS 

Our numerical code is based on the spectral method 
LORENE library routines developed by the Meudon rela- 
tivity group [49j . Our computational grid depends on the 
orbital separation, the mass of the neutron star, and the 
mass ratio. For example, the grid is divided into 10 (8) 
domains for the black hole (neutron star) at large sepa- 

0.15 and mass 
5. Each domain for the black hole 



rations for a neutron star mass of Mg ^ 



ratio MBH/mNS^^^ 

is covered by Nr x Ng x = 41 x 33 x 32 collocation 
points and those for the neutron star by 25 x 17 x 16, 
where Nr, Ng, and N^ denote the collocation points in 
the radial, polar, and azimuthal directions, respectively. 
We use a larger number of collocation points for the black 
hole domains than the neutron star domains because the 
source terms of the gravitational field equations for the 
black hole vary slightly around the neutron star position. 
The amplitude of this excess in the source terms is small 
but significantly affects the accuracy of the results. In 
order to resolve this source in the black hole domains, 
which are centered on the black hole, we need a higher 
angular resolution. On the other hand, the source terms 
of the gravitational field equations for the neutron star 
have large contributions only near the neutron star, and 
drop off monotonically away from the neutron star. Thus 
a smaller angular resolution is sufhcient in the neutron 
star domains, which are centered on the neutron star. 
(See Appendix |X] for the decomposition of the gravita- 
tional field equations into the black hole and neutron 
star parts.) 

As stated above, we choose F = 2 for the adiabatic 
index of the polytropic equation of state in this pa- 
per, and compute several different constant-mass inspi- 
ral sequences. We consider neutron stars of rest mass 




FIG. 1: Contours of the lapse function a in the equatorial 
plane for our innermost configuration of a binary of mass ra- 
tio Mi^^^/Madm.o = 3 and neutron star mass Mf^ = 0.15. 
The minimum value of the lapse, a ~ 0.453, is located on the 
excised surface and the maximum value, a ~ 0.909, at the 
right edge of the figure. The value of the lapse at the point of 
maximum baryon rest-mass density is about 0.558. The con- 
tour curves are spaced by Sa ~ 0.03. The cross "x" indicates 
the position of the rotation axis. 

In Fig.[l]we present contours of the lapse function a for 
a binary of mass ratio M^^/AI^^^^ ~ ^ neutron 
star mass Mg^ = 0.15 at the smallest binary separa- 
tion for which our code converged. The thick solid circle 
on the left-hand side marks the position of the excised 
surface (the apparent horizon), while that on the right- 
hand side marks the position of the neutron star surface. 
The value of the lapse function a on the excised surface 
is not uniform because we use the Neumann boundary 
condition ([M)) rather than a Dirichlet boundary condi- 



6 



tion. The approximate value of the lapse on the excision 
surface is a\s ^ 0.45 - 0.48. 



of 



-0.004 



-0.008 



-0.012 



-0.016 



-0.02, 



BH ,,, NS , 

M. /M.„,.„ =1 
UT ADM.O 




H---nMg =0.10 

3PN approximation 

2PN approximation 

IPN approximation 

Newtonian 



0.02 0.04 0.06 0.08 0.1 0.12 0.14 



-0.003 



-0.006 



0.009 



BH ,, , NS 

/M^„„,„ =3 




• •M|, =0.15 

513'''' = 0.10 

3PN approximation 

2PN approximation 

IPN approximation 

Newtonian 



-0.012 



_0 015 1 • > • > • >^ > • > -r'--!-- l 

0.02 0.04 0.06 0.08 0.1 0.12 0.14 



FIG. 4: Same as Fig. [2] but for sequences of mass ratio 



FIG. 2: Fractional binding energy Eh/Mo as a function of 
QMo for binaries of mass ratio M™/Madm.o = 1- The solid 
line with filled circles shows results for neutron stars of mass 
Mg^ = 0.15, and the dashed line with squares for neutron 
stars of mass of 0.10. The dash-dotted line denotes the re- 
sults of the third post-Newtonian approximation [sol ]. These 
sequences end due to cusp formation - and hence the onset of 
tidal disruption - before the binary reaches the ISCO at the 
turning point of the binding energy. 
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FIG. 3: Same as Fig. [2] but for sequences of mass ratio 
Mi„ /Madm.o = 2. 

In Figs. [2] - [6] we plot the binding energy Eh/Mo 
versus the orbital angular velocity ^IAIq for sequences 
with mass ratios 1, 2, 3, 5 and 10. We find excellent 
agreement with third-order post-Newtonian approxima- 
tions [53, especially for large binary separations. For 
more compact neutron stars with Mg^ = 0.15, we also 
find a turning point in the sequences with mass ratios 

Such turning points mark 
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FIG. 5: Same as Fig. [2] but for sequences of mass ratio 
MS^/MIEm,o = 5. For the more compact neutron star the 
binary now encounters an ISCO before the neutron star is 
tidally disrupted. 



the ISCO, inside of which no stable circular orbits can 
exist. For irrotational binaries this instability is dynam- 
ical (see, e.g. [13]). Binaries with larger neutron star 
compaction and larger mass ratios therefore encounter 
an ISCO before being tidally disrupted. For all other bi- 
naries, the sequences end due to cusp formation - and 
hence at the onset of tidal disruption - before the binary 
encounters an ISCO. 

Qualitatively, this behavior can be understood very 
easily from a simple Newtonian analysis. The binary 
separation at the onset of tidal disruption, dud, can be 
estimated by equating the gravitational and tidal forces 
on a test particle on the surface of the neutron star 
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(43) 
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Same as Fig. [2] but for sequences of mass ratio 

NS 

ADM,0 



= 10. 



If this is smaller than the binary separation at the ISCO, 
at disco /-^BH = 6 for a test particle in Schwarzschild 
coordinates, the binary encounters the ISCO before be- 
ing tidally disrupted. This happens for large mass ratios 
Mbh/Mj^S and large compactions A/ns/-Rns ~ in accor- 
dance with our findings. 

In the test-particle limit the location of the turning 
point should coincide with that of test particle orbit- 
ing a Schwarzschild black hole of mass Mq, ^IMq = 
5-3/2 ^ 0.068. We find similar but slightly larger val- 
ues, riMo - 0.082 and ~ 0.084 for M^^/M^^^^^^ = 5 
and 10. Locating the minimum of a numerically gener- 
ated curve always introduces some additional error, and 
it is possible that this deviation is entirely a numerical 
artifact. However, it is also possible that the deviation 
is a consequence of our "leading-order approximation" of 
the black hole spin. This approximation leads to a small 
but non-zero black hole spin, so that effectively we locate 
the ISCO around a Kerr instead of a Schwarzschild black 
hole. Further evidence for this hypothesis is the fact that 
we do not find turning points of the angular momentum 
that coincide with those of the binding energy. In [Hj, 
the authors similarly found that the binding energy and 
angular momentum for "leading-order" irrotational black 
hole binaries did not have simultaneous turning points. 
In [40], however, the authors show that an improved con- 
dition for nonspinning black holes leads to a very good 
agreement between the binding energy and angular mo- 
mentum turning points. We plan to implement this im- 
proved condition in the near future. 

As one measure of the error in our numerical calcula- 
tions we compute the virial expression (j42p . As shown 
in Fig. [71 all our configurations have relative virial er- 
rors of less than 10"'^, generally on the order of several 
xl0~^. The errors are smaller at intermediate binary 
separations, and somewhat larger both at larger binary 
separations (where the angular resolution becomes worse) 
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FIG. 7: Virial error 5M versus f2A/o for our binary sequences. 



and smaller binary separations (where the larger tidal 
deformation of the neutron star causes larger numerical 
error) . 

From the virial error we can estimate the error in the 
binding energy as follows. We assume that the relative 
error in the ADM mass is of the same order as the relative 
virial error, 



-^ADM — -WaDMj 



.true 



51VL 



(44) 



where AfADM.truo is the true value of the ADM mass. 
From the definition of the binding energy Eq. (|40l) . we 
can write the relative difference of the binding energy 
from its true value i?b.truo as 



Eu — E] 



b,truc 



Eb, 



6M 



,1 



IT 



Mo 



(45) 



Since the term |MADM,truo/(-MADM,true - Mo)\ on the 
right-hand side is of order of 10^ for the binary sepa- 
rations that we calculate, we conclude that the relative 
error in the binding energy is approximately lO^SM . For 
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our 



uui least accurate sequences this results in a relative er- 
ror of several x 10~^ and for our most accurate sequences 
it is a few xlO^'^. 
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TABLE I: Estimated orbital angular velocity at the tidal dis- 
ruption point. 



IlfBH /» ,/-NS 



1 

2 
3 
5 

10 



NS 



QMo 
0.10 M, 



0.0179 
0.0255 
0.0334 
0.0488 
0.0825 



rNS 



0.15 



0.0382 
0.0550 
0.0728 
0.103 



Using a spectral code we cannot follow our sequences 
all the way to cusp formation (see our discussion in [28j). 



To monitor the formation of a cusp we therefore define a 
mass-shedding indicator [4T| 



X 



(a(lnfe)/ar)eq 
(5(ln/i)/9r)poio' 



(46) 



the ratio of the radial derivative of the enthalpy in the 
equatorial plane at the surface to that to the polar di- 
rection at the surface. Given our choice of boundary 
conditions the minimum point of x is not on the X axis 
(which connects the neutron star with the black hole), 
but slightly away from the axis. In Fig. [8] we therefore 
graph the minimum value of x, searched on the surface 
of the neutron star, as a function of SIMq. For spheri- 
cal stars at infinite separation we have Xmin = 1, while 
Xmin = indicates the formation of a cusp and hence 
tidal breakup. Our spectral code stops converging be- 
fore reaching Xmin = 0, but the sharp drop in Xmin is 
an indication of the formation of a cups. Extrapolat- 
ing to Xmin — from the last three points (indicated by 
the thin dotted lines in Fig. [S]) we estimate the orbital 
angular velocity at the onset of tidal disruption (see Ta- 
ble [l|. We did not include the sequence with mass ratio 
M^^/M^^^ f^ = 10 and neutron star mass M^^ = 0.15 
in this analysis, since the binary encounters an ISCO long 
before the neutron star is tidally disrupted (see also Ta- 
ble[lV|. For the sequence with Mj^/^/M^g^.o = 10 and 
Afg^ — 0.10 the extrapolation of Xmin suggests that tidal 
disruption occurs at SIMq ~ 0.0825 (the last entry in Ta- 
ble m , a value very similar to the values of the ISCO in 
those two sequences for which we found a turning point. 
This suggests that the neutron star in this binary may 
be tidally disrupted just as the orbit becomes unstable. 

In '29'|, where we computed similar BHNS binary se- 
quences but with a Kerr-Schild background, we also con- 
sidered a binary of mass ratio M^^/AI^^ = 5 and neu- 
tron star mass = 0.10. There we estimated that 
tidal disruption would set in at ^IMq ~ 0.046. We com- 
pute a sequence with the same definition of the mass ratio 
M^^/M^^ = 5 for comparison, and find tidal disruption 
at approximately flMo ~ 0.053, a value within 15% of 
our earlier value. 

As the orbital separation decreases, the maximum den- 
sity of the neutron star decreases. We define the dimen- 
sionless density parameter 



P 



and monitor the decrease in the maximum density 

_ ^max '?max,0 



^max.O 



(47) 



(48) 



where Qmax.o denotes that of an isolated neutron star with 
the same baryon rest mass. For Afg^ = 0.15 (0.10) we 
have (jniax,o = 0.12665 (0.058827). Note that (5(jmax is 
negative for all configurations considered in this study. 
In Fig. [5] we therefore graph — (Jf/max as a function of the 
orbital angular velocity HMq on a logarithmic scale. 
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FIG. 9: Decrease in the maximum density parameter 5qn 
as a function of Q,Mo. 



Within Newtonian gravity, the decrease in the central 
density scales, to leading order, with 



Sq cx d-" oc 



(49) 



where d is the orbital separation, and for irrotational bi- 
naries the index a takes the value a = 6 (see, e.g. [sij 
for an analytic derivation and [4l[ for comparison with 
numerical results). Here the leading order term is caused 
by a quadrupole deformation in response to the compan- 
ion's tidal field. We see from Fig. [9] that the power-law 
index a is indeed close to 6 for large binary separations. 
For intermediate binary separations, however, the index 
is generally smaller, cr ~ 2 — 3. For the equal-mass bi- 
nary a remains close to 6 for the entire sequence, but 
for increasing mass ratios a decreases at intermediate bi- 
nary separations. We plan to investigate this behavior 
with the help of semi-analytic models in a future publi- 
cation. Finally, for small binary separations, just before 
tidal disruption sets in, a increases again. This has also 
been seen for BNS systems [l^, and may be caused 
by higher order multipole deformations. 
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FIG. 10; The energy spectrum dE^/df for the sequences with 
/M^%yi Q = 1 as a, function of the dimensionless gravita- 
tional wave frequency /Mq. The solid, dashed, dash-dotted, 
and dotted curves show the fits to the Mg^ = 0.15 and 
M^^ = 0.1 sequences and the 3PN and Newtonian expres- 
sions, respectively. Asterisks denote the approximate point 
where a cusp forms and tidal disruption begins, terminating 
the sequence, at frequencies given by Table U 
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FIG. 11: Same as Fig. 1101 but for the sequences with 



Finally, we turn our attention to the energy spec- 
trum of gravitational waves emitted from a BHNS bi- 
nary, which is calculated from the Fourier transform of 
the gravitational wave strains averaged over all angles 
(see, e.g. [12]), 



dE_ 



{Anr')f(\K{fr + \hM)\'' 



(50) 
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FIG. 12: Same as Fig. 1101 but for the sequences with 
M^^ /M^^M O = 5- Note that the sequence with M^^ = 0.15 
reaches a minimum binding energy, where dE^^/df = 0, prior 
to tidal disruption. 



where the brackets denote angle averaging and 



(51) 



where / is the gravitational wave frequency, rather than 
the orbital frequency. From the results of our quasiequi- 
librium sequences, we can compute approximate gravita- 
tional wave energy spectra using the methods described 
in [5^. To do so, we assume that gravitational waves 
from the binary are emitted coherently with f = fl/ir 
at any given moment in time, where f2 is the standard 
orbital angular velocity, and further assume that the sys- 
tem is adiabatic and infall velocities negligible, so that 
the gravitational wave energy losses drive the system 
gradually through a series of quasiequilibrium configu- 
rations. Under these assumptions, we recover a simple 
expression for the gravitational wave energy spectrum. 



dE 
W 



dn 



(52) 



QE 



i.e., we can calculate the energy spectrum by numerically 
differentiating the binding energy along the quasiequi- 
librium sequence with respect to the angular frequency. 
This method has been shown to reproduce closely the ac- 
tual energy spectrum calculated from a fully dynamical 
binary neutron star merger calculation in the frequency 
regime of interest [53| . While we can numerically extrap- 
olate up to the frequency at which we find a minimum 
in the binding energy, representing the ISCO, our results 
do not apply inside the separation at which a cusp forms 
and tidal disruption begins (see Table |T| . Dynamically, 
this would lead to the beginning of mass transfer and 
a rapid evolution away from quasiequilibrium, which in 



turn would produce a qualitatively different gravitational 
wave signal in both the time and frequency domains. 

Instead of taking derivatives of our numerical data for 
the binding energy directly (see Tables HVl and W} . we 
first fit these data to a fourth-order polynomial of the 
form 



Mo 



a^x (1 + aix + a2x'^ + a^x^ + 042;*) (53) 



where x = {^Mop^ = {-KfMof/^. Given this 

parametrization of the polynomial, the coefficients can 

be interpreted as post-Newtonian terms of order i. The 
leading order, Newtonian term, is given by 



(54) 



Here v is the symmetric mass ratio, i.e., the ratio of the 
reduced mass of the system to the total mass. We tab- 
ulate the coefficients a^, together with the exact 3PN 
coefficients for point-mass systems as listed in [50], for 
the different binary sequences in Table [Til We also in- 
clude the "cutoff frequencies" /„ , denoting the frequency 
at which the spectral energy is n% of its Newtonian coun- 
terpart. 



1 dE\, 



ifn) - 



1 dEh 



Mi 



Too 3 



-1/3 



(55) 



for n = 75, 50 and 25. 

In Figs. [To] - [121 we show the results of the energy 
spectrum of gravitational waves. As is expected, bina- 
ries of comparable mass terminate with the formation 
of a cusp, rather than at an ISCO denoting a minimum 
in the binding energy of the binary. For the smallest 
black hole mass we consider, M^^/M^^^^ ~ shown 
in Fig. llOi tidal disruption occurs before we see signifi- 
cant deviations from either the lowest-order Newtonian 
point-mass results or the exact 3PN point-mass expres- 
sion. As the black hole mass and tidal disruption fre- 
quency both increase, so do the deviations in the gravi- 
tational wave energy spectrum away from the Newtonian 
point-mass result. This is especially true for the more 
compact neutron star models with Mg^ = 0.15. For a 
mass ratio M^^/M^^^ ^^ = 5, shown in Fig.m the less 
compact neutron star model with Mg^ = 0.1 (dashed 
curve) shows significant deviations from the 3PN result 
at the tidal disruption point, whereas the more compact 
model with M^^ — 0.15 (solid curve) does reach the 
ISCO before the onset of tidal disruption. This is also 
the case for the most extreme mass ratio we consider, 
^^hr^/^^ADM ~ which any physically reasonably 

compact neutron star model reaches an ISCO prior to 
the tidal disruption point, as had been confirmed by our 
previous dynamical calculations Q. 

We note that the gravitational wave spectrum should 
produce an important constraint on the neutron star 
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compaction, regardless of how the sequence terminates. 
As the compaction decreases, tidal effects grow in impor- 
tance, and tend to increases the binding energy of the 
system (lowering the magnitude of the negative binding 
energy), especially at smaller separations, implying that 
the gravitational wave spectrum will be cutoff at lower 
frequencies. 



TABLE II: Exact 3PN data from 150|], along with polyno- 
mial fit results from the sequences with neutron star masses 
M^^ = 0.15 and 0.10. Coefficients are defined by Eq. 
Characteristic gravitational wave cutofi' frequencies /75M0, 
fboMo, and /25M0 are defined by Eq. ((55}. Cutoff frequencies 
shown in boldface are interpolated from our data, while those 
in regular text are extrapolated. Values shown in italics occur 
at separations within our best estimate of the tidal disruption 
point, given by the values in Table HI and likely overestimate 
the true cutoff frequency. 
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IV. COMPARISON WITH PREVIOUS RESULTS 



We can compare our numerical results with those from 
previous attempts to model BHNS binaries. In partic- 
ular, we compare with our previous work that assumed 
extreme mass ratios [2^ , our comparable-mass results for 
a Kerr-Schild background [2^ , and the results of [lOl for 
comparable masses and a spatially flat background, as in 
this study. 
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FIG. 13: Comparison of our new results to those calculated as- 
■ suming an extreme mass ratio for a spatially flat background 
: geometry (28i] for the minimum of the mass-shedding indica- 
tor Xmin. Here Mbh denotes either the black hole irreducible 
mass (for our new results) or the background mass (for the 
results of jl^); see text. 
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A. Comparison with extreme mass ratio results 



In Fig. [ini we compare our results for a mass ratio 



M^^/M^^ = 10 with those for A/bh/Mb^ = 10 in [2 
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■ for a spatially flat background, but under the assump- 
tion of an extreme mass ratio. Since the tidal effects 
of the neutron star on the black hole are neglected in 

. the extreme mass-ratio approximation, it is impossible to 
.evaluate the black hole's irreducible mass. In [28|, there- 
fore, we instead fixed the background mass of the black 
hole Mbh- Since the difference between and A/bh 
scales with the binary's binding energy (see [sH] for an 

■ analytic, leading-order treatment), the relative difference 

■ is at most a few percent for this mass ratio. 

Our present results agree well with those calculated 
assuming an extreme mass ratio. The difference in the 

■ orbital angular velocity for the same Xmin is order 10%. 
We can explain this 10% discrepancy as follows: for a 
given value of Xmim the tidal deformation of the neutron 
star by the black hole is approximately the same in the 
two cases. This implies that the orbital separation is very 
similar as well. However, the orbital angular velocity in 
our calculation here is given in terms of the distance from 
the center of mass of the binary system to the center of 
the neutron star, while for an extreme mass ratio it is 
given by the binary separation. This is because with the 
extreme mass-ratio assumption, we assumed that the bi- 
nary's center of mass is located at the center of the black 
hole. The difference between these two definitions of dis- 
tance, along with a change in the black hole mass we 
use in the definition of the mass ratio, accounts for the 
discrepancy in the orbital angular velocities. To Newto- 
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nian order, there is a difference between = \J Mq / in 
our calculation liere and = ^Mbh/^"^ for the extreme 
mass ratio. The relative difference between these angular 
velocities is about 10% for the mass ratio in question. 



TABLE III; Comparison between our implementation and 
that of GOG. Here "Ex. BC" stands for the excision boundary 
condition, and "Dec. of ((T} - ([9])" stands for the decomposi- 
tion of the relevant equations and variables. 



B. Comparison with Kerr-Schild coordinates 
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FIG. 14; Comparison of our new results to those calculated 
using a Kerr-Schild background ,29] for the total angular mo- 
mentum. The dash-dotted line denotes the results of the third 
post-Newtonian approximation [sol ]. 

In Fig. 1141 we compare our new results for the total 
angular momentum with those in poj . for a mass ratio 
M^^/M^^ = 5 and neutron star mass M^^ = 0.10. The 
difference between our treatment here and that in j2^ is 
the choice of the spatial background metric jij . In this 
pap er, we assume a flat background, jij = r]ij, while in 
[29j we chose a Kerr-Schild metric. We see from Fig. [14] 
that the behavior of the sequence agrees to within 5% 
[5^ throughout most of the sequence. 



C. Comparison with previous flat background 
results 

In [s^l (hereafter G06), Grandclement computed sim- 
ilar sequences of BHNS binaries, adopting an approach 
that in some ways is similar to ours. Both codes are 
based on the Lorene spectral method library, but are 
otherwise completely independent. In particular, we im- 
pose different inner boundary conditions on the black 
hole's excision surface, use different conditions to make 
the black hole nonspinning, and adopt a different decom- 
position of the equations and variables. We summarize 
these differences between our implementations in Table 

imi 
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FIG. 15; Comparison between our results and those of G06 
for the binding energy as a function of QMo, for a binary 
of mass ratio M^^/M^^yi g = 5. The solid line with filled 
circles marks our results for a neutron star compaction of 
Af^DMo/^o = 0.150, which corresponds to the baryon rest 
mass of Mg^ — 0.153. The dashed line with triangles rep- 
resents those of G06 for a neutron star compaction of 0.150. 
The dash-dotted line denotes the results of the third post- 
Newtonian approximation [50l |. 



To compare our results with those of G06, we show in 
Fig. [15] the binding energy as a function of QMq for two 
sequences with the same physical parameters. Both se- 
quences are for a binary with mass ratio M^^/M^^j^ g = 
5. Our results are for a neutron star rest mass of 
Mg^ = 0.153, which corresponds to a compaction of 
-^^ADM 0/^0 = 0.150. We find agreement with the third- 
order post-Newtonian result to within 3% except for con- 
figurations inside the turning point, and hence inside the 
ISCO. We also include the binding energy for a neutron 
star compaction of 0.150 as computed by G06 (and as 
found in Fig. 3 of GOG). Clearly this curve shows a 
much larger deviation from both our results and the post- 
Newtonian ones, by about a factor of two. We speculate 
that this deviation could be caused by the larger numer- 
ical error in the calculations of G06. According to G06, 
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the relative virial error is as large as 2%. If the masses 
themselves carry an error of 2%, the binding energy may 
have an error as large as 100% or more, as we have out- 
lined in the paragraph including Eqs. (I44p and ([^5]). This 
is consistent with the deviations found in Fig. ]T5\ 

V. SUMMARY 

We construct quasiequilibrium sequences of BHNS bi- 
naries in general relativity. We solve the constraint equa- 
tions of general relativity, decomposed in the conformal 
thin-sandwich formalism and subject to the black hole 
boundary conditions of [1^ , together with the relativistic 
equations of hydrostatic equilibrium. In contrast to our 
earlier approach in (29| , where we solved these equations 
for a Kerr-Schild background, we now adopt a spatially 
flat background. We also employ a new decomposition 
of the equations and variables that leads to a significant 
improvement in the numerical accuracy. 

We construct constant-mass sequences of mass ratios 
M^^/Mf^^Q = 1, 2, 3, 5, and 10 for neutron star 

baryon rest masses of Mq^ = 0.15 and 0.10. We find 
excellent agreement between our calculated binding ener- 
gies and the third-order post-Newtonian approximation, 
especially for large binary separations, but also agreeing 
to within 5% at small separations. We locate the onset 
of tidal disruption and, for large mass ratios and neutron 
star compactions, the ISCO. 

We find some evidence that our "leading-order approx- 
imation" for constructing nonspinning black holes leads 
to small errors, and plan to implement an improved con- 
dition in the near future. 
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APPENDIX A: DECOMPOSITION OF THE 
GRAVITATIONAL FIELD VARIABLES 

For the construction of binary configurations with 
codes that are based on spherical polar coordinates it 
is natural to use two computational domains, each one 
centered on one of the binary companions. An equation 
of the form 

Act) = R-HS (Al) 
can then be split into two equations 

A0BH = RHSbh (A2) 



A</.NS = RHSns, (A3) 

where (j) = ^bh + '/'ns and RHS = RHSbh + RHSns- The 
two equations (jA2p and (|A3P can then be solved on two 
computational domains, one centered on the black hole 
and the other centered on the neutron star. 

Clearly, the separation of the source terms RHS into 
its two parts RHSbh and RHSns is not unique. One 
guiding principle is to move those parts of RHS that are 
large in the neighborhood of the black hole into RHSbh , 
and likewise for the neutron star. Another principle is 
that each of the source terms should asymptotically co- 
incide with that for the isolated black hole and neutron 
star when the orbital separation is large. In addition, 
we use this decomposition freedom to deal with the fact 
that metric quantities are not defined inside the excised 
parts of the black hole. It is not clear how to treat such 
functions if they appear in RHSns • One approach is to 
artificially fill the resulting "hole" in the function with 
a smooth function (see, e.g. |57[). Instead, we carefully 
separate the source terms in such a way that no metric 
quantities that are excised in the black hole interior affect 
the right-hand side of the neutron star equation. 

Specifically, we first decompose the metric quantities 

as 

a^aBH+am, (A4) 
P'^pBH+(3m+Kou (A5) 

^Z" = '0BH + ^NS, (A6) 

where (31^^. is the rotational shift vector. We then decom- 
pose the gravitational field equations as follows. Equa- 
tion ([7]) becomes 

+l{^m + c^)-'A^fA^4s^ (A7) 
A^NS = -27r^V- i(V^NS + c^r' AfA^s^ {A8) 
Eq. ([5]) is decomposed as 

APbb + ^5'(a,/3^H) = 2i'^a,(a^-6) 

-2A^sdj[{am + Ca)(^Ns + c^)"^], (A9) 

+2i^ga,[(aNS + c„)(i/'NS + c^)"^], (AlO) 
and Eq. ^ is written as 

AckBH = a^p-^A,JA'^ - 2r]'^diadj IniP 

-("NS + Ca)(V'NS + C4,y^Afj^A^g 

+2t]'^ di{aj^s + Ca)dj ln(V'NS + c^), (All) 
AaNS — 47ra'0^(p + 5*) 

+(aNS + Ca)(V'NS + c^y^AffA% 
-2-q'^di{ans + Ca)9j ln(V'NS + c^,), (A12) 
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where we define A^^g as 

^ ^ (^NS + c,)^ /^.^^ ^ ^ _ 2 .,5^,^. A . 

(A13) 

Note that the neutron star equations (|A8p , (jAlOp , and 
(|A12|) contain metric quantities Hke a and that are not 
defined inside the excised black hole, but only in prod- 
ucts with matter quantities that vanish outside the neu- 
tron star. As a consequence, the neutron star equations 
are regular everywhere. All singular, or un-defined terms 
have been moved to the black hole equations, where they 
are harmless since they are excised from the computa- 
tional grid. 

We finally point out that the outer boundary con- 
ditions of the total lapse function and total conformal 
factor are asymptotically flat, i.e., a\r^oc — 1 and 
^p\r^oc = 1- This implies that each part of the metric 
quantities, i.e., aBH, Q;nSi ^bh, and ^/^ns, cannot go to 



unity individually at infinity. We set the outer bound- 
ary conditions for these quantities as aBnlr^oo = 0.5, 
ckNslr^oo = 0.5, V'BHlr^oo = 0.5, and V'Nslr^oo = 0.5 for 
the convenience of the computation. We therefore insert 
constants Ca — 0.5 and = 0.5 to ensure that the total 
quantities take on their proper asymptotic values. 



APPENDIX B: TABLES OF SEQUENCES 

We summarize our results in Tables ITVl and fVl In these 
tables, we tabulate the orbital angular velocity 17, bind- 
ing energy iJb, total angular momentum J, decrease in 
the maximum density parameter (Jgmax, minimum of the 
mass-shedding indicator Xmin, and fractional difference 
SM between the ADM mass and the Komar mass along 
a sequence. 
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TABLE IV: Physical parameters for a binary sequence with neutron star baryon rest mass M^^ = 0.15. The ADM mass, and 
the isotropic coordinate radius of the neutron star in isolation arc Madm.o = 0.139 and fo — 0.815 (k = 1). The neutron star 
compaction is M^dm o/^o = 0.145 where Rq is the areal radius of the isolated neutron star. We list the binding energy E^,, 
total angular momentum J, orbital angular velocity Q, decrease of the maximum density parameter Sqmax, minimum of the 
mass-shedding indicator Xmin, and fractional difference SM between the ADM mass and the Komar mass. The prefix f denotes 
the position of the turning point. 



Mass ratio: Af?^/MADM,o = 1 

d/Mo QMo Eb/Mo J/M'o ggmax Xmin SM 

15.28 0.0153 -7.34(-3) 1.107 -9.56(-5) 0.961 6.66(-6) 

13.10 0.0191 -8.36(-3) 1.048 -4.16(-4) 0.934 2.83(-5) 

10.91 0.0246 -9.71(-3) 0.985 -1.45(-3) 0.871 1.95(-5) 
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8.51 0.0347 -1.18(-2) 0.914 -7.80(-3) 0.575 4.76(-5) 

8.30 0.0359 -1.20(-2) 0.907 -1.03(-2) 0.435 1.42(-4) 

Mass ratio: M?"/Madm,o = 2 

d/Mp nMp E^/Mq J/Mj dq^^ Xmin SM 

20.37 0.0101 -5.15(-3) 1.089 -3.65(-5) 0.992 4.45(-4) 

11.64 0.0224 -8.30(-3) 0.891 -1.19(-3) 0.952 1.33(-4) 
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TABLE V: Same as Table IIVI but for the neutron-star baryon rest mass M^^ = 0.10. The ADM mass, and the isotropic 
coordinate radius of the neutron star in isolation are M^dm.o = 0.0956 and fo — 0.990 (ft = 1). The neutron star compaction 
is M^g^,o/-Ro = 0.0879. 
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